Many-body localized molecular orbital approach to molecular transport 
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An ab initio based theoretical approach to describe nonequilibrium many-body effects in molecular 
transport is developed. We introduce a basis of localized molecular orbit als and formulate the 
many-body model in this basis. In particular, the Hubbard- Anderson Hamiltonian is derived for 
single-molecule junctions with intermediate coupling to the leads. As an example we consider a 
benzenedithiol junction with gold electrodes. An effective few-level model is obtained, from which 
spectral and transport properties are computed and analyzed. Electron-electron interaction crucially 
affects transport and induces multiscale Coulomb blockade at low biases. At large bias, transport 
through asymmetrically coupled molecular edge states results in the occurrence of "anomalous" 
conductance features, i.e., of peaks with unexpectedly large/small height or even not located at the 
expected resonance energies. 

PACS numbers: 73.63.-b, 85.65.+h 
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I. INTRODUCTION 

Experimental and theoretical investigations of elec- 
tron transport through single molecules are in the fo- 
cus of the rapidly growing field of molecular electron- 
ics^— . Electron-electron interaction plays an important 
role, controlling the position of resonant levels and lead- 
ing to phenomena such as Coulomb blockade and Kondo 
effects. Depending on the ratio between the energy scales 
associated with an effective charging energy and coupling 
to the leads, molecular junctions can be classified in three 
groups. In the case of very small coupling, the molecular 
orbitals are only weakly disturbed, strong charge quan- 
tization and Coulomb blockade take place and transport 
is mainly determined by sequential tunneling^ - —. In the 
opposite case of large coupling, the electronic molecular 
orbitals are hybridized with states in the leads, charge 
quantization is suppressed, transport is mainly coherent 
and the conductance is of the order of the conductance 
quantum Go = e 2 //*£— . Finally, in the intermediate sit- 
uation, the effective electronic spectrum of a molecule is 
determined essentially by the hybridization, and the in- 
terplay between charge quantization and coherent trans- 
port may be important^ — . Let us consider as a com- 
monly used benchmark example a gold-benzenedithiol- 
gold (Au-BDT-Au) molecular junction with the central 
molecule S-CqH^-S. The experiments^^ show that, 
while it is difficult or impossible to observe the typi- 
cal Coulomb blockade features in this type of molecular 
junctions, there are signatures of correlated transport, 
namely a conductance gap at small voltages and a com- 
plex shape of the conductance peaks. One can conclude 
that transport in the case of intermediate coupling can 
be correlated and partially incoherent. 

In parallel with the experimental investigations, a 
number of theoretical methods were suggested to de- 



scribe the structure and electronic properties of molec- 
ular junctions. In the case of coherent transport ab ini- 
tio DFT+NGF methods^"—, combining density func- 
tional theory (DFT) and nonequilibrium Green function 
(NGF) techniques^— , have become standard^. How- 
ever, the use of DFT, which is a powerful tool to deal with 
ground state properties, may become questionable when 
applied to transport and nonequilibrium situations, espe- 
cially when inelastic and interaction effects are essential. 
Indeed, the use of the solutions of the Kohn-Sham equa- 
tion as physical quasiparticle states can not be rigorously 
established. Besides, the DFT+NGF is a mean-field the- 
ory and lacks to describe many-body effects in systems 
with strong electron interactions. Thus, being based on 
the basis of an atomistic modeling of molecular junctions, 
the rigorous extension of the DFT approach to describe 
transport phenomena remains a challenge. 

On the other hand, model-based approaches are par- 
ticularly important to understand the physics of corre- 
lated transport. The models are solved usually within 
two main approaches: the NGF technique in the case 
of strong coupling to the leads, and the quantum mas- 
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Figure 1: (Color online) Schematic picture of a Au-BDT-Au 
molecular junction. The dashed box defines the extended 
molecule. It comprises the central molecular region and parts 
of the leads, also marked with boxes. 



2 



ter equation (QME) 25 i 26 in the weak-coupling limit. The 
quantum master equation is usually formulated in the 
basis of many-body eigenstates of the molecule. It gives 
a fairly complete description of sequential tunneling, the 
main features of Coulomb blockade and even can cap- 
ture Kondo physics for temperatures of the order of or 
larger than the Kondo temperature^ 7 .. Finally, for some 
nonperturbative effects covering the whole range of tem- 
peratures, e.g. the Kondo effect in the crossover regime, 
more sophisticated treatments are necessary, e.g. numer- 
ical renormalization group approache s 28 ! 29 , other numer- 
ical method s 30 ! 31 or Keldysh field theories 3 - 2 -. 

The challenge for molecular transport theory is to com- 
bine an ab initio approach, required to take into ac- 
count a realistic geometry and capable to provide the 
electronic structure of molecular junctions, with a many- 
body quantum transport technique, which is necessary 
to incorporate the correlation effects. In this direction a 
number of different ab initio based approaches were sug- 
gested 3 - 3 ^—. It should be noted that many-body calcula- 
tions usually require sophisticated analytical and numer- 
ical methods and can be very time consuming. Hence the 
methods of transport theory can not usually be applied 
directly to large realistic systems; instead a combined 
approach is preferable, where an effective model takes 
into account only the states predominantly participat- 
ing in transport. Progress in this direction was achieved 
recently in applications to Coulomb blockade phenom- 
enal, correlated transport through atomic wire s 35 i 36 , 
many-body interference 37 and Kondo effect^— . Be- 
sides, Coulomb blockade phenomena in benzene and ben- 
zenedithiol junctions were also considered on the basis 
of semiempirical atomistic models^—. However, a sys- 
tematic ab initio based many-body theory of molecular 
transport is still lacking. 

Our aim is to further develop such a theory. The 
main problem to be solved on the way from an atom- 
istic model to transport calculations is that a huge num- 
ber of microscopic (single-atom) degrees of freedom must 
be reduced to an effective few-electron (or few many- 
body state) interacting model, as a prerequisite of suc- 
cessive many-body transport methods. The main build- 
ing blocks of our approach are an effective electron model 
(an orthonormal basis of localized molecular orbitals and 
a many-body Hamiltonian in this basis) extracted from 
atomistic calculations, and a nonequilibrium quantum 
transport method, based on nonequilibrium Green func- 
tions or on the quantum master equation. 

In the case of strong or intermediate coupling to the 
leads the electronic molecular levels should be considered 
together with the levels in the leads. A corresponding 
generic atomistic model is shown in Fig.[TJ A so-called 
extended molecule (inside the dashed box) is placed be- 
tween equilibrium electrodes. The extended molecule 
consists of the central region (roughly the molecule) and 
two leads (the regions of electrodes near the molecule). 
The size and boundaries of the central region are actually 
not fixed and should be determined in a way to include 



all relevant electronic states, as we will see below. 

After having defined the appropriate size of the ex- 
tended molecule, we proceed as follows: 

i) We perform Hartree-Fock or DFT ab initio calcula- 
tions within the extended molecule and obtain the molec- 
ular orbitals, which are orthogonal and can serve as a 
basis for a many-body Hamiltonian. 

ii) The central region (e.g. an organic molecule) and 
the metallic leads have quite different physical properties, 
and the approximations required to describe interactions 
are also different. Thus, it is convenient to transform the 
basis of molecular orbitals obtained in (i) into a basis of 
localized molecular orbitals (LMOs), which can be spa- 
tially separated into a basis for the central region and a 
basis in the lead ends. Besides, the advantage of LMOs is 
that the Coulomb interaction in this basis can be simpli- 
fied to the Hubbard form. This procedure enables us to 
obtain a many-body Hubbard- Anderson Hamiltonian in 
the central region with ab initio calculated parameters, 
including the coupling to the leads. The rest of the leads 
can be treated within some mean-field approximation. 

iii) Using the Hubbard- Anderson Hamiltonian, many- 
body methods (nonequilibrium Green functions and the 
master equation in this paper) can be effectively applied 
to analyze spectral and transport properties of molecular 
junctions. 

The paper is organized as follows. The ab initio ba- 
sis of localized molecular orbitals is introduced in Sec.HIl 
The electron-electron interaction and Hubbard approxi- 
mation are discussed in Sec.lIIIl The parameters of the ef- 
fective Hubbard- Anderson Hamiltonian for the subspace 
of electronic states relevant for transport through the 
central region are derived in Sec. HVI In Sec.|Vl we use 
nonequilibrium Green functions to analyze the coupling 
to the electrodes and spectral properties. The many- 
body spectrum of the central region is discussed in Sec. [VI] 
and the transport results are shown in Sec. lVIII We give 
a short conclusion in Sec. I Villi 



II. LOCALIZED MOLECULAR ORBITALS 

The first stage of our approach includes ab initio calcu- 
lations of the optimized geometry and the relevant basis 
of electronic states. For the preliminary geometry opti- 
mization and molecular dynamics of whole structures we 
apply the DFT code Siesta^. For the extended molecule 
only (without the full electrodes) the full-electron quan- 
tum chemistry code Firefly (former PC GAMESS)^ is 
used. The final geometry optimization is performed us- 
ing a hybrid DFT method, usually B3LYP. Then, the 
calculation of the molecular orbitals (MOs) is performed 
inside the extended molecule (including both the central 
region and leads). Our test calculations show that a sim- 
ple local density approximation (LDA) gives reasonable 
results for the considered systems. At this stage we ob- 
tain the MOs ^ MO (r) with associated energies Ei. MOs 
have the advantage of being normalized and orthogonal. 
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Figure 2: (Color online) The energy spectrum of the molec- 
ular orbitals of the extended molecule shown in Fig. [I] The 
occupied levels are marked red, the empty levels yellow. 



However, the canonical MOs of the extended molecule 
can not serve as a good basis for systems with interac- 
tions, because both electron-electron and electron- vibron 
interactions require physically different approximations 
in different, metal or molecular, parts of the junction. 
It is hence more convenient to use localized molecular 
orbitals (LMOs). The interactions between localized or- 
bitals have a simple and transparent form and the appro- 
priate approximations can be found. Besides, the num- 
ber of required, relevant basis states is smaller and better 
controlled for LMOs than for MOs. 

To proceed, we divide all MOs into four groups (Fig. [2]). 
Most relevant for transport are the transport levels near 
the Fermi energy of the electrodes. These levels are se- 
lected for the localization procedure and include both oc- 
cupied and valence MOs in an appropriate energy range 
around the Fermi energy. This energy range should be 
larger than the energy scales of the external bias volt- 
age and temperature. The other criterion is that the 
obtained LMOs should be localized strongly enough in- 
side the central region and in the leads: only in this case 
it is possible to separate the system into parts and use 
different approximations for the central region and leads 
separately. Alternatively, the partial localization of only 
relevant MOs (for example only 7r-type orbitals) is possi- 
ble. In any case, these transport electron states play the 
main role. All other polarization states, which are fur- 
ther away from the Fermi energy or do not participate 




Figure 3: (Color online) Localized molecular orbitals in the 
central region of an Au-BDT-Au molecular junction. 



in transport for some other reasons, can still affect the 
interaction between transport electrons and, in partic- 
ular, screen the Coulomb interaction. The polarization 
MOs can be localized in the same way as the transport 
orbitals. Finally, the core orbitals at low energies can be 
included into the effective (pseudo)potential and empty 
orbitals at high energies are neglected. 

For benzenedithiol, considered in this paper, the trans- 
port window was chosen to be about 4 eV, and contains 
about 60 MOs, mainly due to the large density of states 
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Figure 4: (Color online) Localized molecular orbitals in the 
leads (several out of many are shown by different colors). 

in the metal leads. The parameters of the obtained ef- 
fective model are rather robust against the exact choice 
of this number. 

The LMOs are obtained from the MOs by the unitary 
transformation 

V^°(r) = Y>^ MO (r). (1) 

i 

The indices with bars a,/3... denote the states without 
the spin degree of freedom. We apply the Foster-Boys 
localization method 53 , which minimizes the spatial ex- 
tent of the orbitals and maximizes the distance between 
orbital centers. Thus, we obtain maximally localized or- 
bitals. Out of the about 60 orbitals only 5 are localized in 
the central region. While the initial MOs spread across 
the extended molecule, the LMOs are spatially localized 
in the central region (Fig. [3]) or in the leads (Fig. [4]). 

Due to the unitary transformation (pQ) the LMOs are 
still orthogonal and normalized, but the expression 

i 

is no longer diagonal. Moreover, the related Hamiltonian 
H LMO takes the form 

/ H L Vlc \ 

jjLMO = xLMO , (3 ) 

V V RC Hr J 

where Hl, Hq MO , and Hr are the Hamiltonians of the 
left lead, the central region, and the right lead separately. 
The direct coupling between left and right leads can be 
neglected in most cases, because the LMOs of different 
leads are only very weakly overlapping. 

III. COULOMB INTERACTION AND 
HUBBARD APPROXIMATION 

Having the LMOs at hand we can calculate the 
Coulomb matrix elements. The electron-electron inter- 
actions are described by the Hamiltonian 

^ ee = 2 V ^,isdid^d s d 7 , (4) 



where a = {a,a a } and a a denotes the spin. The matrix 
elements for the scalar Coulomb interaction V(|r — r'|) 
are defined as 

Vaj3,j5 = %,7^<7 a a 7 ^tr65 (5) 

V & - M s = J dvj dv'r & {v)r B {v')V{\v-v'\)^{v)^{v'), 

(6) 

where S aa(J/3 is the Kronecker symbol. For the systems 
with localized wave functions V ; a( r ), where the overlap 
between two different states is weak, the main matrix 
elements are those with a = 7 and (3 = 5. We checked 
that, indeed, the overlap of 3 or 4 different orbitals can 
be neglected in comparison with the overlap of only 2 
orbitals. In these cases it suffices to replace H ee by the 
Hubbard Hamiltonian 

H H = \ ^ U a ph a hp, (7) 

describing only density-density interactions with 
n a — d^d a and the Hubbard parameters defined as 

U aP = J dvj dr'|^ MO (r)| 2 |^°(r')| 2 n|r-r'|). 

(8) 

As a further advantage of LMOs, the local nature of elec- 
tron correlation is better described in the localized basis. 

The interaction V(|r — r f \) is the bare Coulomb inter- 
action V(r) = 1/r, if we include all electrons into the 
localization procedure. Actually we restrict the effec- 
tive Hamiltonian to transport electronic states, on which 
we performed the Hubbard approximation. The remain- 
ing polarization electrons are included only through a 
screened Coulomb interaction. The screening can be 
described by the effective interaction potential in RPA 
(or GW) approximation, which is, however, energy de- 
pendent. To keep the simplicity of the Hubbard ap- 
proximation, we use a screened Coulomb interaction 
V(r) = l/(er) with e w 1.5. This approximation gives 
reasonable values of U a p for the benchmark 7r-orbital 
model of benzene CeH6 which we compared with the opti- 
mized semi-empirical PPP model^. In this way all coeffi- 
cients are derived, but further semi-empirical corrections 
could be included for better agreement with experiments. 



IV. MANY-BODY MODEL 

The next important step is the derivation of the en- 
tire effective Hamiltonian in the basis of the LMOs for 
the central region. As we already explained, we divide 
the extended molecule into the central part (actually the 
molecule in our particular case) and the leads (Fig. [5]). 
The full Hamiltonian is the sum of the noninteracting 
molecular Hamiltonian H^, the electron-electron inter- 
action Hamiltonian i^ ee , the Hamiltonians of the ends 
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of the leads Hr(l), an d the tunneling Hamiltonian 
describing the molecule-to-lead coupling: 



H = H% + H ee + Hl + Hr + #r- 



(9) 



In our case the central Hamiltonian Hq + # ee is replaced 
by the Hubbard cluster Hamiltonian: 

Hc = H% + H H = J2t*pdid(3 + \ y ^U OL ^h OL n^ (10) 



where e a p = e a /3 + e(p a ^a(3 are the bare energies of the 
LMOs, including the shifts due to an external voltage. 

The noninteracting molecular Hamiltonian is ob- 
tained from the LMO Hamiltonian H^ MO , Eq. (j3]). The 
zero- voltage energies € a p = €^S aa(Tp are calculated from 
the HF or DFT MO energies E^ from which the contri- 
bution of the Hartree terms due to Hubbard interactions 
should be substracted: 



6 ap = S od EiS Ae aaS t 



(ii) 



where 



= J2J2 m 2 \ u ** s *i + 2U ^ 1 - • ( 12 ) 

i 7 

In this expression n% denote the populations of the cor- 
responding LMOs in the ab initio calculation. The sum 
is taken over all occupied molecular orbitals. 

The coupling to the leads is described by the tunneling 
Hamiltonian 

H T = Y.{ V sk<rAkJa + V ska , a dic ak ,) , (13) 

s=L,R ka,a 

and the Hamiltonians of the left and right leads are 

H S =L(R) = ^2e ska cl ka c ska: (14) 

ka 



Left lead 



Central region 



Right lead 



H, 



GO H c Q0- H R 




Figure 5: The model of the extended molecule. 



where k is the index of a state, and a is the spin. Note 
that in our case the states in the leads are not plane 
waves, but are represented by LMOs, calculated simulta- 
neously with the LMOs in the active region. The leads 
are considered at the mean-field (DFT) level. The equi- 
librium electrodes, which can have different electrochemi- 
cal potentials, determine the boundary conditions for the 
leads. 

For the 5-level model (Fig. [3]), which represents actu- 
ally 10 levels with spin, we obtain the following parame- 
ters (for spin up or down, all numbers are in eV): 

/ -17.80 0.06 0.00 0.03 0.06 \ 
0.06 -19.81 -0.04 -0.54 -0.02 
{e^} = 0.00 -0.04 -18.79 0.04 -0.06 
0.03 -0.54 0.04 -18.20 0.00 
\ 0.06 -0.02 -0.06 0.00 -16.25 / 

(15) 

The matrix of the Hubbard parameters calculated from 
expression ([8]) reads 



{UaB} 



/ 4.08 2.72 2.31 1.36 1.36 \ 

2.72 4.32 2.43 2.34 1.36 

2.31 2.43 3.84 1.59 1.59 . (16) 

1.36 2.34 1.59 4.00 2.72 

V 1.36 1.36 1.59 2.72 4.08 / 



The diagonal entries in this matrix correspond to the in- 
teraction of two electrons in the same LMO state but 
with different spin. The off-diagonal terms denote cou- 
pling beween two different LMOs spin of the electronic 
state. Of course, these expressions should be transformed 
into the 10-level basis before performing further calcula- 
tions. 

At finite bias voltage V (defined by the left and 
right electrical potentials, V = ipL — cpr) the energies 
are shifted. In linear approximation these shifts are de- 
scribed by r] a : e al3 = e al3 + e(p a S a ^ <p a = <p R + rj a ((p L - 
<Pr), where the parameters < rj a < 1 characterize 
the symmetry of the voltage drop across the junction, 
rj a = 0.5 stands for the symmetric case. Note that the 
energies and other parameters can also be ab initio calcu- 
lated at finite voltage, but that is very time-consuming. 

The coupling matrix elements V ska ^ a in Eq. ([T3|) are ob- 
tained directly from the localization procedure. Indeed, 
the Hamiltonian of an extended molecule takes the form 
Eq. (j3]), and the off-diagonal terms describe the coupling 
to the leads. The number of states in the leads is many 
times larger than in the central region. Thus, to leading 
approximation, we can average over the lead level dis- 
tributions and couplings (so-called wide-band limit). In 
this approximation the level- width function 



r s>Q! /3(e) = 2ir^V sk(T ^V* k(J ,J(t - e ska ) 



(17) 



is energy independent. 

We now return to our 5-level model. The couplings of 
the first two states (localized at the left sulfur atom) and 
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the last two states (localized at the right sulfur atom) 
are characterized by the level- width functions T^ii = 
0.16 eV, r L2 2 = 0.21 eV, T RM = 0.28 eV, T R55 = 0.1 eV. 
All other couplings are small and do not play an es- 
sential role. Thus, all parameters of the model Hamil- 
tonian Eqs. (fT0|) - (fT4|) are well defined and we can pro- 
ceed with transport calculations. In view of the exper- 
iments [lO. 11], we will perform calculations below room 
temperature, /c^T < 0.025 eV, implying T > ksT. 

As we discussed in the introduction, transport at finite 
voltage can be described in the framework of nonequilib- 
rium Green function or quantum master equation ap- 
proaches implying numerical methods. For benzene- 
based junctions several methods were used, including 
coherent DFT based, the master equation approach in 
the sequential tunneling limi t 33 i 48 , sophisticated approx- 
imations in the framework of the nonequilibrium Green 
function metho d 40 i 46 i 47 i 49 , and other methods 3 ^. In this 
paper we use both NGF and QME methods, trying to 
attack the problem from both sides. It should be noted, 
however, that for our case, ksT <T <C U, both NGF and 
a QME with second order rates can only give a qualita- 
tive description of the transport problem. Very recently, 
a QME approach for a single-level junction able to de- 
scribe the regime T ~ has been proposed 55 . An 
extension of this method to a multilevel system will be 
subject of future investigations. 



V. NONEQUILIBRIUM GREEN FUNCTION 
APPROACH TO SPECTRAL PROPERTIES 

We follow the formulation pioneered by Meir, 
Wingreen and Jauho — The lesser (retarded, ad- 
vanced) Green function matrix G < ( jR ' A ) = G^^'^ of 
a nonequilibrium molecule can be found from the Dyson- 
Keldysh equations in the integral form 

G R (e) = G*(e) + G${e)± R (e)G R (e), (18) 

G<(e) = G R (e)t<(e)G A (e), (19) 

or from the corresponding equations in the differential 
form 



A 

7/3" 



(21) 



Here 



£R,A,< = £R,A,<T) + £R,A,<T) + £R,A,<(I) (22) 

is the total self-energy of the molecule composed of the 
interaction self-energy E jR ' A ' < ( / ) and the tunneling (cou- 



pling to the left (L) and right (R) lead) self-energies 

yR,<(T) _ T >f /yR,< 
^s=L,R 

,fl,<(T) _ ^ J G Ri< V. 



^-Tr J =V] c G^<V sC , (23) 



J sotf3 



ka 



{ V s*ka,a G fka Vaka,?} , (24) 



where Gf^ ,<: is the Green function of the leads. 

The retarded tunneling self-energy E^ T ^ can be rep- 
resented as 

T \e) = A s (e - e Vs ) - l -f s (e - e Vs ), (25) 

where A s is the real part of the self-energy, which usually 
can be included in the single-particle Hamiltonian H^, 
and T s describes level broadening due to coupling to the 
leads. In the case of noninter acting leads with continuous 
energy spectra, the level- width function is determined by 
the expression ([TT|) . 

For the corresponding lesser function of the noninter- 
acting leads one finds 

£< T \e) = it s (e - e Vs )f s {e - op.), (26) 



where 



exp((e-/x s )/T) + l 



(27) 



is the equilibrium Fermi distribution function with chem- 
ical potential ji s (ks = 1). 

The expression for the interaction self-energy can not 
be obtained exactly. In the nonequilibrium Hartree-Fock 
approximation one has 59 



v R(i) 



(n 7 ) 5 a p, 



(28) 
(29) 



We do not consider more sophisticated cases here. 

The current from the left (s = L) or right (s = R) lead 
into the central system is described by the expression (for 
spin-unpolarized leads) 



s=L,R 



^Tr{f s (6-e^) (G<(6)+ 



f°(e-eip s )[G R (e)-G A (e)])} 



(30) 



Applying the NGF technique to our case, we should 
take into account that our system initially consists not 
of three, but of five parts: the large electrodes, quantum 
leads and the central region (Fig. [I]). Accordingly, the 
full Hamiltonian has the form 



H 



Vt 1 ' h l V lc 



V 



el] 

Vlc 






t 

RC 



o \ 






V RC H r vf 



(31) 
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function of the central region, 



Figure 6: (Color online) The spectral function within the dif- 
ferent approximations: initial DFT (black dashed), restricted 
HF (red), unrestricted HF (green) and equation- of- mot ion 
(dashed blue). 



where the central part is the same as before, Eq. (J3j), and 
the additional terms describe the electrodes and the cou- 
pling between the electrodes and the leads. 

The solution of Eq. (jTHJ) for the central part is in this 



fiR 



(e + irj)I ■ 



H° - 



(32) 



with the lead self-energies X 
width functions defined as 

n {T) = n C G R L VLc, 



R(T) 



(see Eq. (j23j)) and level- 



Y T = -2Im£ 



R(T) 



Y R = -2Im£ 



R(T) 



(33) 
(34) 



The lead Green functions for noninteracting leads (or for 
leads described by the effective mean-field Hamiltonians 
H s ) are defined correspondingly by the electrode self- 
energies 



1 



(e + iri)i-H L -tf 
1 

{e + i n )I-H R -tf 



(35) 
(36) 



The calculation of the electrode self-energies is done by 
standard methods^. The lesser Green functions are cal- 
culated in the same way. We assume additionally that 
the distribution function in the leads is the same as in 
the corresponding electrodes. 

The equations are solved self-consistently within four 
approximations: initial DFT with the mean- field ener- 
gies Eq. (j2j), the restricted HF approximation (RHF) with 
n&t — n ot^ the unrestricted HF approximation ([23]) and, 
finally, the equation of motion method (EOM ) 59 i 6Q . 

First we analyze the equilibrium (zero bias) spectral 



A(e) = -2^ImG 



R 

Coca. ' 



(37) 



The first thing one can see (Fig. [6]) is that the RHF ap- 
proximation gives a spectral function similar to the ab 
initio (DFT) one. This is not surprising as the ab ini- 
tio calculation is also RHF and we simply extracted the 
Hartree contribution when calculating the energies e^p 
in Eq. (fT2j) . The other important point is that the re- 
sults obtained in HF and EOM approximations are dis- 
tinctly different and a gap is opened at the Fermi surface. 
The analysis of the populations n a = (h a ) of the single- 
particle states shows that two empty states are located 
at the left and one at the right side of the molecule (sec- 
ond and fifth states in Fig. [3]). These two empty states 
have the same spin (in the HF approximation the ground 
state is spin polarized and degenerate, but quantum fluc- 
tuations can switch between different spin orientations), 
indicating that the true ground state, with quantum fluc- 
tuations taken into account, can be spin singlet or triplet. 
We discuss this point in the next section. 



VI. THE MANY-BODY SPECTRUM: 
EXACT DIAGONALIZATION AND THE 
GROUND STATE PROPERTIES 

To gain insight into the many-body energy spectrum of 
the central system we perform an exact diagonalization 
of Eq. ([TO]) obtaining the set of many-body eigenstates 
| A). Calculating the tunneling matrix elements we obtain 
from Eqs. ([T0|) and ([T3]) the Hamiltonian 

H c + Ht = Y,Ex\\)(\\ (38) 

A 

ska + T* kaX , x 



s=L, J R;fecr:AA / 



^^W*<A'l4|A>. (39) 



First, we analyze the many-body spectrum. With 5 
LMOs we get 1024 many-body states in the Fock space. 
The lowest 8-particle states consist of a series of alternat- 
ing singlets and triplets (see Table|I]). In particular, the 
ground state is a singlet, practically degenerate with a 
triplet {E& — E% g w 10 _4 eV). It follows, at a distance 
of roughly 0.3eF a second pair of singlet-triplet quasi 
degenerate states. The 9-particle states are all doublets 
with a relatively regular distance of the order of 0.5eV. 
Finally, it is important to keep in mind that the energies 
of the lowest four 8-particle levels lie all below the one 
of the 9-particle ground state. The 7-particle states have 
much higher energies. 

Note that in the sequential-tunneling master equation 
method the exact many-body states can be partially oc- 
cupied at finite temperatures, but not at zero tempera- 
ture, and the level broadening is not taken into account. 
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Level 


Energy [eVJ 


ibpm [n\ 




n 1 1 a n 
-91.1049 


U 


8' 


-91.1848 


1 


8" 


-90.8653 





8'" 


-90.8648 


1 


% 


-90.7866 


1/2 


9' 


-90.3693 


1/2 


9" 


-90.0891 


1/2 



Table I: The eigenenergies and the associated spins of the 
lowest 8 and 9 particle levels of an Au-BDT-Au molecular 
junction. 



This can give some noticeable difference in the position 
of the transport resonances compared to the NGF calcu- 
lations, where the levels can be partially occupied even at 
low temperatures, and where the real part of the HF self- 
energy, Eq. ([23]) , describes the energy shift of the single- 
particle levels. 

The population probabilities P\ are found from the 
master equation 



(40) 



where the tunneling rates are, in second order in the tun- 
neling Hamiltonian, 



^xx' 



]T [l&.fiiEx-Ex.-etp.) 

s=L,R]<j 

+7>7x (1 " fi{Ex> -E X - op.))] , (41) 



with 



Ixx 



-ry^T sk(Ti \\'T* k(TXX ,5(E\ - E\> - e ska 
a k 

_ EW(A|4|AKa(A'|4|» -E X , -e 



2tt 



ska J 



a(3k 



Y r -,^ (E X - E\r ) (A 1 4 1 A') ( A' I d a | A) . (42) 

a(3 



This expression connects the tunneling rates to the level- 
width function; thus the T S(Tja ^(e) calculated by the NGF 
method can be used, see Eqs. ([33134ft . In the wide-band 
limit one has hj xx , = 27rpo|^scr,AA / 1 2 > where po is the 
density of states. 

To check that the simple (diagonal) form of the master 
equation can be used, we have analyzed the many-body 
spectrum of the considered system and came to the con- 
clusion that no coherences are needed for the description 
of the transport since the degeneracies are not of orbital 
but of spin nature (e.g. triplets for the 8-particle and dou- 
blets for the 9-particle states). However, there cannot be 
mixing of states with different total spin since otherwise 
the mixing will depend on the choice of the direction of 
the quantization axis. The solution of the Eqs. ()40ti42|) is 



straightforward and can be obtained by direct numerical 
integrations in stationary and time-dependent cases. 

As we discussed before, in the equilibrium state at zero 
voltage there are 8 electrons distributed due to thermal 
smearing between the states 8^ and 8', see Tab. HI An 
equilibrium occupation with 8 electrons is in agreement 
with the HF calculations of Fig. [6] In Tablellfl the com- 
position of the many-body states in terms of the five lo- 
calized molecular orbitals of Fig. [3] is quantified in terms 
of the average populations n a = (h a ) of the single par- 
ticle states obtained from exact diagonalization of He 
(see Eq. ([TO]) ). The composition of the states 8^ and 8' 
is similar to the HF average populations. As discussed 
in Sec. lVIII and shown in Fig.[8j the LMO occupations 
change at finite bias. 





ni 


n 2 


n 3 


714 


n 5 


NGF 


1.9725 


1.0934 


1.9989 


1.9294 


1.0220 


8, 


1.9801 


1.1337 


1.9973 


1.8685 


1.0204 


8' 


1.9798 


1.1339 


1.9972 


1.8685 


1.0207 


8" 


1.0390 


1.8047 


1.9997 


1.1753 


1.9815 


8 /// 


1.0269 


1.8159 


1.9995 


1.1764 


1.9812 


9, 


1.9990 


1.4657 


1.9992 


1.5365 


1.9996 


9' 


1.9510 


2.0000 


1.9989 


1.9999 


1.0502 


9" 


1.0772 


1.9917 


1.9998 


1.9780 


1.9532 



Table II: The average populations n& of the single-particle 
states at zero bias voltage. Calculations within the NGF 
method are shown in the second line. They agree with the 
composition of the states 8 g and 8 X obtained from exact diag- 
onalization of He (see Eq. {TO])). 



VII. TRANSPORT CHARACTERISTICS 

We are now able to calculate and interpret the current- 
voltage characteristics of the benzenedithiol junction. 
The current at finite voltage, which is given by 

Is=l,r = e bTx'f°(Ex ~ E\i - ev s )P x , 

AA';<r 

- 7 & (1 - f a {Ex ~ E\i - op.)) Px] , (43) 

is presented together with the differential conductance 
in Fig. [7] The curves are asymmetric with respect to a 
bias inversion because the junction geometry was cho- 
sen to be slightly asymmetric. As the main result we 
find a multi-scale Coulomb blockade. The large region of 
suppressed current is about 2 Volt wide. However, the 
current is completely blocked only in a much smaller re- 
gion of bias voltage, as small steps in the current (peaks 
in the conductance) are present at lower biases. 

As a first step in the analysis of the current voltage 
characteristics we consider the average particle number in 
the central system presented in the left panel of Fig. [8] At 
low biases the average particle number is 8 correspond- 
ing to the neutral configuration of benzenedithiol. The 
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V[V] 

Figure 7: (Color online) Current-voltage curve (black solid 
line) and differential conductance (red dashed line). 



many-body state with the minimal grand-canonical en- 
ergy (Eg = E—fiN) is in fact the 8-particles ground state 
(see Fig. [TO]). When the bias drop is raised in the junction 
the average particle number takes values between 8 and 
9 ensuring that the dominant transitions are negative ion 
resonances. 

A further insight into the dynamics is obtained by 
monitoring the average occupation of the different lo- 
calized molecular orbitals shown in the right panel of 
Fig. [8] At low biases the symmetric central orbital (the 
third orbital from the top in Fig. [3]) is completely oc- 
cupied, ns = 2. Its occupation undergoes a sensible 
variation only at the voltages of the large current steps 
Vb « ±1.5V. Large variations in the population of the 
asymmetric LMOs centered around the molecule-lead in- 
terface (orbitals 1,2,4 and 5 in Fig. 3) are instead associ- 
ated to the small current steps present at lower voltages. 
Interestingly, at a bias of Vb « 0.5 V the effective spatial 
symmetry of the system is recovered with the populations 
of the asymmetric states being all equal. 

A deeper understanding of the dynamics of the sys- 
tem is obtained by the analysis of the occupation of the 
many-body states (Fig. [9]), their energies (Tab. [I]), and 
the transition rates among them (Tab. lIII( schematically 
represented in Fig.fTO]). If the calculation of the current 
is performed taking into account hundreds of many-body 
states, the essential physics at the biases presented in 
Fig.[7]is captured by considering the lowest four 8-particle 
levels (for a total of 8 states) and the lowest three 9- 
particle ones (6 states). 

The tunnelling events from (to) the source or the drain 
connect these many-body states. The tunnelling rate 
T AA is the product of a geometrical part (73^/ of Eq. (|32J) ) 
and an energetic contribution encoding the energy con- 
servation in the tunnelling event and the Pauli exclusion 
principle (see Eq. flSJ). The energetic contribution en- 
sures that the rate T xx changes (and correspondingly 
the current through the system) every time the resonant 



Particle number LMO populations 




-2 2-2 2 

Figure 8: Bias dependence of the average electron number 
(left) and average individual populations (right) on the molec- 
ular junction. 

condition E\ — E\> — e<j) s = is fulfilled. With this argu- 
ment it is already possible to assign a specific transition 
to most of the peaks in the conductance of Fig. [71 In par- 
ticular the transitions 8 P , 8' O 9' are associated with the 
peak at the most negative bias and 8^ o 9^ to the second 
peak from the left. The first small peak at positive bias 
is anomalous and we will return to it later. We only note 
that its position depends on the temperature and that 
it moves to the 8^ O 9 g resonance at low temperatures. 
The rightmost conductance peaks are instead associated 
to the transitions 8", 8"' o 9' and 8", 8'" O 9", respec- 
tively. 

The approximate symmetries of molecular geometry 
are very important since they introduce selection rules 
which distinguish between transitions which are energet- 
ically equally allowed. In Tab. lIIII we report the transition 
rates 7 scr between the different many-body states. Here 
the values are given in eV and the spin a is chosen to 
fulfill spin conservation in the tunnelling event. In the 
case of a doublet to triplet transition the value of the 
rate reported is the one involving the triplet state with 
maximum projection along the quantization axis. Except 
for the transition 8' «->• 9" all transitions show a very pro- 
nounced left-right asymmetry. It is much easier for ex- 
ample for an electron to tunnel in (or out) of the molecule 
from (to) the left instead of the right lead when this tun- 
nelling event involves the many-body eigenstates 8' and 
9'. This asymmetry is essential to explain the dynamics 



N = 8 N = 9 




Figure 9: The occupation of the manybody states 
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Figure 10: The energies, tunneling rates and the associated populations of the most relevant states (see discussion in the text). 



1 L ° 


% 


9' 


9" 




0.0017 


0.1025 


0.0024 


8' 


0.0056 


0.2039 


0.004 


8" 


0.1033 


0.0003 


0.0074 


8 „, 


0.1676 


0.0013 


0.0488 






9' 


9" 




0.0442 


0.0127 





8' 


0.0866 


0.0278 


0.0022 


8" 


0.0044 


0.0056 


0.1157 


g /// 


0.0013 


0.011 


0.2185 



Table III: The transition rates 7 s 
manybody states. 



between the different 



sures than even after the opening of the current channel 
the occupation of the 8^ (together with the almost de- 
generate 8') is still the largest one. In the right panel 
of Fig. [10] we schematically represent the tunnelling rates 
and the associated populations of the most relevant levels 
for a bias just above (in absolute value) the first negative 
bias conductance peak. Starting from this population 
distribution it is then natural to observe the next visible 
current step related to the transition 8^(8') o 9'. Since 
this time the left tunnelling rate dominates, the popu- 
lation of the 8-particle states decreases substantially in 
favor of the 9-particle ones. Generally, a more uniform 
mixing of states with different particle number is associ- 
ated with a larger fluctuation of the number of electrons 
in the central system and thus with a larger current. 



of the system at low biases and can be understood in 
terms of the spatial distribution of the many-body eigen- 
states. 

The left transition rate is larger than the right one 
when the transition from an 8 to a 9-particle state is as- 
sociated to a larger variation of the density in the orbitals 
1 or 2 than in the orbitals 4 or 5. Analogous arguments 
hold for the reverse situation. 

Let us now return to the interpretation of the current 
voltage characteristics with the help of Fig. [TO] By con- 
vention, to a positive bias voltage Vb corresponds a sta- 
tionary particle current flowing from right to left while 
the electrical current flows from left to right. We concen- 
trate first on the negative bias. From an accurate analysis 
of the definition of the tunnelling rates (Eq. flH}) it is not 
difficult to prove that the first step in the current is due to 
the resonant condition between the S g (8 f ) and 9 g states 
at the left lead. Current flows since the system oscillates 
between the 8 g (8 f ) and 9 g states by receiving an electron 
from the left lead and by releasing it to the right one. The 
asymmetry between the transition rates, j Ra > 7 Lcr , en- 



The dynamics at positive bias is more complex. In 
particular the first conductance peak occurs at a bias at 
which even the ground to ground state transition is not 
yet open. This anomalous behaviour is understandable 
when taking into account the large left-right asymme- 
try of the rates. As schematically represented in the left 
panel of Fig.[10j even before the (right lead) resonance 
between the 8 g (8') and the 9 g state opens a conventional 
current channel, the states 8" and 8'" get strongly pop- 
ulated. The fundamental reason is the large probability 
to tunnel out of the system at the left lead through the 
transition 9 g — >> 8", 8'" which is also energetically favor- 
able. Very soon the states 8" and 8'" become the new 
effective ground states for the system (see Fig. [9]). In 
this scheme it is thus not surprising that i) the first con- 
ductance peak is located at an "average" between the 
8 g {8') 9 g resonance and the 8" (8'") O 9^ one; ii) the 
next two conductance peaks at positive bias occur at the 
8"(8'") <H> 9' and 8" (8'") <r> 9" resonant conditions. 
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VIII. CONCLUSIONS 

In conclusion, we developed a many-body localized 
molecular orbital approach to transport through molec- 
ular junctions with the following protocol: 

1. Geometry optimization using DFT and hybrid 
DFT (usually B3LYP based) methods. 

2. Molecular vibrons can be calculated after the geom- 
etry optimization (not considered in this paper). 

3. Molecular orbit als of the extended molecule are ob- 
tained. Localized molecular orbitals (LMOs) are 
constructed and form the basis for all subsequent 
calculations. 

4. A Hubbard interaction is introduced for the LMOs 
in the central region: only density- density Coulomb 
integrals are taken into account. 

5. Electron- vibron interaction can be included in the 
central region (not considered in this paper). 

6. The leads are kept as effectively noninteract- 
ing (mean-field approximation). The interac- 
tion Hamiltonian between leads and central region 
yields the relevant tunneling couplings. 

7. A spectral analysis and transport calculations are 
performed on the basis of the ab initio based 
Hubbard- Anderson model. 

Using the benchmark example of a benzene-dithiol 
molecular junction, we performed the full line of calcula- 
tions in the framework of this approach. We determined 
the geometry of the junction, calculated molecular or- 
bitals and transformed them into localized molecular or- 
bitals. Upon using an energy range of about 4 eV around 
the Fermi energy of gold, we obtained a basis of 5 LMOs 
with energies e a p. Then we calculated the Coulomb ma- 
trix elements U a p for these orbitals and coupling matrix 
elements V s k a ,a between the central region and the leads. 
Using the parameters e a ^, U a p and V s ko-,a, obtained from 
ab initio calculations, we calculated the spectral function 
in the framework of the nonequilibrium Green function 
approach (in the RHF, HF and NEOM approximations). 
Besides, the model was transformed into the many-body 



eigenstate basis, and the quantum master equation (ap- 
plied in the sequential tunneling limit) was used to cal- 
culate the current. It is shown that transport through 
asymmetrically-coupled molecular edge states results in 
suppressed peaks of the differential conductance at small 
voltage and unexpectedly large peaks at higher voltages. 
The origin of these anomalies could be explained upon 
analyzing the occupation probabilities of the many-body 
states as well as their compositions in terms of LMOs. 
In general, we could qualitatively understand the equi- 
librium state and main transport properties of the con- 
sidered molecular junction with strong electron-electron 
interaction and intermediate coupling to the leads. 

Nevertheless, the further development of the theory is 
necessary with respect to both ab initio and quantum 
transport aspects. The results presented in this paper 
are only partially self-consistent because the parameters 
tap, U a f3 and V s k a ,a are calculated at zero voltage, but 
used at all voltages. Actually, it is possible to extend 
the theory to include the recalculation of the parameters 
at finite voltage and the influence of the nonequilibrium 
charge in the central region on the leads. A related issue 
is the effect of the external field on the LMOs energies, 
which we treat using a simplified linear approximation. 
The Hubbard interaction plays the main role, but the 
corrections due to non density-density interactions and 
polarization of the molecule can be important as well. 
Finally, we expect that the method proposed in Ref. [HH 
could be of importance to treat the parameter regime 
/c#T < T < U typical for molecular junctions with inter- 
mediate coupling to the leads. 
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